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Abstract 

An accreting massive black hole (MBH) in a galactic nucleus is surrounded by a dense stellar cluster We 
analyze and simulate numerically the evolution of a thin accretion disk due to its internal viscous torques, due 
to the frame-dragging torques of a spinning MBH (the Bardeen-Petterson effect) and due to the orbit-averaged 
gravitational torques by the stars (Resonant Relaxation). We show that the evolution of the MBH mass accretion 
rate, the MBH spin growth rate, and the covering fraction of the disk relative to the central ionizing continuum 
source, are all strongly coupled to the stochastic fluctuations of the stellar potential via the warps that the stellar 
torques excite in the disk. These lead to fluctuations by factors of up to a few in these quantities over a wide 
range of timescales, with most of the power on timescales > {M,/Mj)P{Rj), where M, and Mj are the masses 
of the MBH and disk, and P is the orbital period at the disk's mass-weighted mean radius Rii. The response of 
the disk is stronger the lighter it is and the more centrally concentrated the stellar cusp. As proof of concept, 
we simulate the evolution of the low-mass maser disk in NGC 4258, and show that its observed (9(10°) warp 
can be driven by the stellar torques. We also show that the frame-dragging of a massive AGN disk couples the 
stochastic stellar torques to the MBH spin and can excite a jitter of a few degrees in its direction relative to that 
of the disk's outer regions. 

Subject headings: galaxies: nuclei — galaxies; active — black hole physics — stars: kinematics and dynamics 
— accretion disks — galaxies: individual (NGC 4258) 



1. INTRODUCTION 

Massive black holes (MBH) gain most of their mass, 
and presumably a substantial fraction of t h eir spin, in the 
course of luminous accretion (lSoltan|[l982t lYu & Tremaind 
12002V This implies accretion via a radiatively efficient, 
geometrically thin and optically thick accretion disk. As 
the MBH mass grows, and in the absence of strong exter- 
nal perturbations (e.g. galactic mergers), the stellar popu- 
lation around it is expected to settle into a centrally con- 
centrated cusp, r ij, with 0.5 < 7 < 2.5, whether adi- 
abatically (e.g. I Yound 119801). or by two-body re laxation 
! (lBahcall & w5l II 977t l Alexander &Hopmanll200^ This 
stellar cluster, which spatially coexists with the disk and 
\ extends farther out to a substan tial frac t ion of the MBH's 
radius of dynamical influence jMerritll 120061 Fig. 12), 
can affect the disk both directly , via hydrodynamical in- 
\ teractions (e.g. 'Syer et at 119911 : lArtymowicz et alJ 119931 : 
IVilkoviskij & Czernyu2002 ), and indirectly, by pu rely grav- 
itational interactions jBregman & Alexander! I20()9l hereafter 
BA09), and thereby possibly influenc e the cosmic evolution 
of MBHs. As we argue below (Section [T3] ). unless the disk is 
very compact, it is the gravitational interactions that exert the 
dominant effect by deforming the disk's geometry. 

The best known example of a deformed thin Keplerian ac- 
cretion disk is that of the circumnuclear maser disk of ac- 
tive galaxy NGC 4258, which displays an (9(10°) warp on 
the 0.1 pc scale. In a previous preliminary test of concept 
study (BA09), we analyzed the physical scales in the prob- 
lem, and argued that the observed warp is consistent with 
resonant torquing of the disk by random fluctuations in the 
stellar distribution around the MBH. Here we extend and gen- 
eralize that study by simulating numerically the evolution of 
a thin accretion disk due to its internal viscous torques, due 
to the frame-dragging torques of a spinning MBH and the 



stellar orbit-averaged gravitational torques. We revisit the 
NGC 4258 maser disk, and further explore the implications 
of stellar torquing for an accretion disk's geometry, mass ac- 
cretion rate and covering factor, and for the spin evolution of 
MBHs. 

The paper is organized as follows. In the remainder of 
the introduction we describe a simplified model for a galac- 
tic nucleus harboring an accreting MBH, and scale it to the 
NGC 4258 system (Section fTTTTi. W e then briefly review Res- 
onant Relaxation (RR, Section [L2b an d de rive some relevant 
order of magnitude estimates (Section [1.31 ). In Section |2] we 
describe how the MBH / disk / stellar cluster system is mod- 
eled and simulated numerically (a brief summary of the nu- 
merical scheme is presented in appendix |A]i. We present the 
results in Section [3] and discuss and summarize them in Sec- 
tion [H 

1.1. The NGC 4258 system 

Maser-emitting nuclear accretion disks are unique, clean 
probes of the environment near MBHs. The first discov- 
ered and best-studied maser disk, NGC 4258, is also the 
thinnest and most Keplerian (to better than 1%; iMalonevI 
120021) nuclear disk yet discovered. Radio observations of 
H2O maser emission from the edge-on disk ( iGreenhill et alj 
11995b reveal the disk morphology, gas velocities and ac- 
celerations, and allow accurate measurements of the MBH 
mass (M, = 3.7 x lO^M©), the distance to the host galaxy 
(D = 7.2 ± 0.3 Mpc) and the sp atial extent of the maser regio n 
{Ra = 0.13 pc to /Jfo = 0.26 pc) jHerrnstein et al.|[7996l [19991) . 
Together with optical and X-ray observations, these also con- 
strain the disk's accretion r ate (10^^ < M < lO^^Moyr"') 
( iNeufeld & M alonev! [T99 5I). The NGC 4258 m aser disk 
shows a clear (9(10°) warp ( jHerrnstein et aljl996l) . It is note- 
worthy that observations also indicate possible warps in the 
maser disks of Circinus ( iGreenhill et al.,2003.) and NGC 3393 
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jKondratko et al.ll2008l) . 

Here the stellar cusp around the MBH is modeled 
as a single mass power-law cusp of stars, A^*(< r) ~ 
Nh{r/rijf'^^, which extends inward from the radius of in- 
fluence at r/, = GMt/af, where the velocity dispersion of 
the galactic bulge O), c an be estimated from the empiri- 
cal M, I a relation ( Ferr arese & Merritlll2000t iGebhardt etaTI 
l2003HShields et al.ll2003k and contains Nh = llhM»/Mi, stars 
of mass M* each, where pih ^ 0(1). We model NGC 4258 
with r/, = 7 pc, and assume jU/, = 2 (the formal value for 
an MBH-less singular isothermal distribution), = IMq, 
and u nless stated otherwise, 7= 1.75 (a dynamically relaxed 
cusp, iBahcall & Woig[T976l) . For the purpose of simulating 
the torques on the maser disk of NGC 4258, only stars in the 
radial range 0.01 pc to r/, are considered. The maser disk 
model is described in detail in Section 1X21 



1.2. Resonant relaxation 



Resonant relaxati o n dRauch & Tremaind 119961 : 
iRauch&IngaiisI [19981 iHopman & Alexanded l2006h is a 
rapid angular momentum relaxation mechanism that operates 
in potentials with a high degree of approximate symmetry, 
which restricts orbital evolution (e.g. fixed ellipses in a nearly 
Keplerian potential, or fixed planar rosettes in a nearly spher- 
ical potential). In such potentials the fixed, orbit-averaged, 
stellar mass distribution exerts a constant residual torque on 
a test mass, which persists over a coherence time fo as long 
as perturbations due to deviations from the perfect symmetry 
remain small. The accumulated change in angular momentum 
J over fo then becomes the "mean free path" in J-space for the 
noncoherent random walk phase on timescales longer than fo- 
When fo is long, the mean free path is large and the random 
walk is rapid. The efficiency of RR is determined by the 
nature of the physical process that perturbs the symmetry and 
limits fo (e.g the Keplerian symmetry is perturbed far from the 
MBH by the potential of the stars and near it by relativistic 
precession). For circular orbits in a near spherical potential, 
such as those of gas streams in an accretion disk, the RR 
torques can only change the orbital orie ntation, but not the 
eccentricity, i.e. J{r) = Jc{r) = const (iGlirkan & HopmanI 
l2007h . where Jc = \/GM,r is the maximal (circular) angular 
momentum at radius r. On timescales f < fo, the direction of 
the angular momentum vector of a circular orbit of radius r 
changes coherently due to the residual forces by stars on the 
same scale, 

w(r) = |A7x(r)| A.(r) - ^^^K{r){M,/M,)t /P{r) , (1) 

where A'* is the number of stars within r and P is the or- 
bital period. The numeric prefactor )3x ~ 0{1) may de- 
pend somewhat on the parameters of the cusp, and can be 
determined from A^-body simulations, where it is measured 
to be jSx > V2 (lEilonet al.ll2009t Kupi & Alexander 2011, 
in prep.); j3x — a/2 is adopted here. The "warp factor" w 
(0 < vv < 2) corresponds to a tilt in J by cos / = 1 — 

In the limit where the disk mass M^i is negligible com- 
pared to the MBH mass M,, the coherence time is set by the 
randomizing effect of RR itself on the torquing stars ("self- 
quenching"). In this case, the coherence time for the torque 
on a circular orbit of radius r by stars on the same scale is 



'sq 



(r)^Asq(M./M,)P(r)/V^^, 



(2) 



of A^-body simulations indicates thatAsq = 1.0 ±0.1 (Kupi & 
Alexander 201 1, in prep.); Asq = 1 is adopted here. The warp 
factor grows over the self-quenching time to n'(fsq) = l5±Asq ~ 

\/2, that is, quenching occurs when the accumulated tilt in the 
orbital direction grows to ; ~ 7t/2. An approximate correction 
for the coherence time when the back reaction of the disk on 
the stars cannot be neglected, is introduced in Section [T3] 



1.3. Order of magnitude estimates 

We argue here, based on an approximate analysis of scales, 
that purely gravitational interactions between the stars and 
the disk, rather than hydrodynamic ones, typically domi- 
nate the torquing of the disk; that around lower mass MBHs 
(< 10^ Mq ) the disk warps in response to the stellar torques, 
while around more massive MBHs it changes its orientation 
as a nearly rigid body; and that stellar torques are expected to 
warp the maser disk of NGC 4258 by 0(10°). 

Direct stars / disk interactions — Stars crossing the disk exert 
a torque on it by direct hydrodynamic interaction. To assess 
whether this effect competes with the purely gravitational RR 
torques, we estimate the torque exerted on the disk in reaction 
to the ram pressure the disk exerts on the stars. The residual 
torque by A^* stars of mass and radius Ri,, on randomly 
oriented orbits, in a volume of size R enclosing a disk of mass 
Md, is 



' RFram ^ RVN^pv llr] 



X J 



(3) 

where p is the disk density, v" ^ GM,/R, rx = ma\{R^,rB) 
is the effective cross-section radius of the star, where rg = 
2GMi,/{v^ + c]) ~ 2{M^/M,)R is the Bondi radius in the 
limit Cs < O(10kms-') < v - O(10^kms-'), which holds 
for a cold accretion disk close to the MBH. For = IMq, 
R^ = lRr.„ and for M. = 3.7 x 10^ M©, R = 0.26 pc and 
H/R = 0.002 (the values for NGC 4258, iHerrnstein et"an 
l2005h . rx — R*. The gas density can be estimated by p '--^ 
Mdl{TtR^2H), where H is the disk's exponential scale height. 
We further assume that the cold thin disk is close to its self- 
gravity stability limit (as implied by the non-s mooth morphol- 
ogy of many obs erved galactic maser disks, iMalonevI 12001 
iBraatz & Gugliucci 2008. and is pre dicted to be the general 
case in AGN disks, lGoodmanll2003l) . M^/M. - H/R. The 
ram torque is then 



^ M,t GM, , I ^GM, 
2kR^H R * 2 R 



while the residual RR torque on the disk by the stars is 



Trr 



GMi^M, 



d 



R 



(4) 



(5) 



The RR torques affect the disk continuously over an orbital 
time, whereas the ram torque is applied only twice per orbit, 
when the star crosses the disk. The ratio between the dura- 
tions the torques are active is Afram/AfRR ~ AH / {2nR), and 
therefore the ratio of accumulated changes in the disk angular 
momentum per orbital time is 



Atn 



A/rr Trr AfRR 



1 M. 



(6) 



where Asq is an (9(1) numeric prefactor. PreUminary analysis 



For the parameters of NGC 4258, AJam/^JRR ^ 10"^. Even 
if AGN disks are limited by gravitational instabiUty to radii as 
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small as - 2000^^ (r^ = GM./c^), as argued by iGoodmaiil 

(l2003b . then A/ram/A/RR - O.O2(M./lO^M0)-\ and RR 
torquing dominates for all MBHs with mass M, > few x 
IO^Mq. The torquing of a thin disk by ram pressure is there- 
fore negligible, and the torquing effects of the stars on the disk 
are well approximated by pure gravitational interactions. 

Stellar torques vs. disk angular momentum transport — The ef- 
fect of the RR torques on the disk geometry depends on the 
ratio between the stellar torques and the rate of angular mo- 
mentum diffusion in the disk. When the internal viscous 
torques are small, differential torquing across the disk re- 
sults in warps, which persist on the warp diffusion timescale 
^warp ~ {R/H)^P /2Ka2, where ^ 0{l) is the dimension- 
less viscosity parameter associated with the vertical viscosity 
(Eq. \T2\ . The necessary condition for the disk to be Keple- 
rian, A^^, <C Mt/M-^, and the condition for the warp to persist 
against diffusion, tQ < fwarp, provide together a soft upper limit 
to the mass range of MBHs where disk warping is efficient. 



(7) 



where A^* = /atM./M* is the number of stars out of the 
Ni, ^ 0{M,/Mi,) in the MBH radius of influence that are close 
enough to the disk to affect it appreciably. Conversely, when 
to > fwarp- the disk flattens out rapidly, and the stellar torquing 
results in an overall evolution in the inclination of the disk 
as a nearly rigid body. The upper limit for warping depends 
quite sensitively on the assumed parameters, especially the 
disk's aspect ratio. For R/H = 250, /at = 0. 1 , = 1 Mq and 
a2 = l,it isM. -C»(1O^M0). 

Residual stellar angular momentum — A random density per- 
turbation on the spatial scale r has a Poisson magnitude 
N-t (r) and carries angular momentum of the order Jn ~ 

yjNi,{r)Mi,^GM,r. This is the maximal angular momen- 
tum that can be transferred to the disk from that scale. The 
torque exerted by the stars on the disk will lead to an equal 
and opposite torque by the disk on the stars. By the time 
the disk has been torqued by Jn, the reaction force on the 
stars will disperse the perturbation, and another, uncorrected 
one, will take its place. The back-reaction thus introduces 
a new coherence timescale for the perturbation, freact('") such 
that |AJ(freact)| = JN{r), which may be shorter than the self- 
quenching coherence timescale, fjq and thus limit the action 
of the torques on the disk. 

We estimate fieact here by approximating the disk, which 
extends between Ra and R^, as a thin ring of mass Mj 

and mass-weighted mean radius Rj ~ InM^^ J^^R^Y-AR, 
and evaluating the stellar torque on the disk in the lim- 
its of small and large spatial scales. The lever arm on 
the gas disk is Rj. When r ^ R^i, the magnitude of the 
force by the stars on the disk is y/Njj)GM-i,/r^, and so 
AJ{t)/JM = ^i_A,^{Rd/ryi'-tlh^{r). Conversely, when 
r ^ Rd, the force on the disk is y/Njj)GMi, / R^^, and so 

M{t)/Jc{Rd) = ^i_A^ci{r/Rd?/^t/ts^{r). In shorthand nota- 
tion, 

AJ{t)/J,{Rd) = j3xA,q(r/7;,/)i/2-0f/^(^) , (8) 
where = sign{r — R^). The back-reaction timescale is then 



(9) 

which increases monotonically with r. The coher- 
ence time is to{r,Rd) = min(4q,freact)- Since typically 

tveact{r,Rd)/h<l{r,Rd) = yK{r}{M,/Md)/V2 < 1 for r ~ R,,, 
where the stellar torques on the disk are most effective, the 
back-reaction time is what sets the limit on the coherence 
time. The change in the angular momentum of the disk over 
the back-reaction time is 



Wreactl'") 




(10) 



For the NGC 4258 cusp parameters (S ection [TTTI ) and a disk 
mass ofMj - 3000 (lMartinll2008l: BA09) with r ^ R^ 
0.16 pc, the predicted overall change in disk orientation is 
~ 27°, and the differential change across the inner and outer 
edges of the masing region, the observed warp angle o, is an 
order unity fraction of that (BA09). We therefore anticipate 
that RR warping can lead to the (9(10°) warp that is observed 
in NGC 4258. 

2. CALCULATIONS 

2.1. The evolution equation 

The equation governing the evolution of the angular mo- 
mentum surface density, L, of a thin accretion disk under the 
influence of internal and external torques, is given in the limits 
of a Keplerian velocity field, no azimu thal modes and diffu- 
sive warp propagation (a\ >H/R ), by (iPapaloizou & Pringl^ 
[1981 iPringlell 1 992t IOgiivia[T999h 



dL 1 d 

'dt'^R'dR 

]_d_ 
^R'dR 

]_d_ 
^R'dR 



V2R- 



dR 



1 dt 
-V2RL— 

2 dR 



V,RL X - 



(11) 



where L ~ L^/GMtR, i ^h/L and E is the disk's mass sur- 
face density. The disk is thus effectively modeled as a set 
of concentric rigid thin annuli. The first term on the right- 
hand side of Eq. ( fTTT i describes the angular momentum that is 
carried to the central sink with the inflowing mass, the sec- 
ond term the angular momentum that is advected outward 
by the disk's internal viscous torques, which are expressed 
by the azimuthal kinematic viscosity Vi (responsible for the 
mass inflow) and the vertical kinematic viscosity V2 (respon- 
sible for the unwarping of the disk). The third term de- 
scribes the precessional to rque that accompanies large ampli- 
tude warps (IOgilvielll999h . We assume isotropic viscosities, 
V„ = OLnCiH (n = 1,2,3), where c, is the mid-plane isothermal 
sound speed, and use the second order expansion of the three 
dimensionless viscosity parameters 



1,2,3) (12) 



in terms of the local dimensionless local warp amplitude 
\^ ^R \dtldR\ (.Ogilvift.1999; .Lodato & Price. .20 10.) . The 
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first order coefficients are 



2(l+7af) 
ai (4 + af ) 



3(1 -2ar) 
2(4 + a2) 



(13) 



We find that the 2nd order terms do not change the results sig- 
nificantly in our simulations. The second order coefficients 
are listed for completeness in appendix |2.3.1l For typical val- 
ues ai < 1, V2 ^ Vi, that is, the disk unwarps faster than it 
flows into the MBH. 

Text = Tbp + Trr expresses the external torques per unit 
area applied to the disk, which here are those due to relativis- 
tic frame-dragging and to the stellar torques. Tsrc expresses 
the torque per unit area due to a source term at the disk's 
outer edge (Section |2. 3. lb . The effect of frame dragging on 
the disk, the Bardeen-Petterson (BP) effect, is given to lowest 
post-Newtonian order b}{3 (iBardeen & Pettersonll 1975b 

TBP = flLTXL, (14) 

where Qlt = {2G/c^R^)J, is the Lense-Thirring precession 
angular frequency and J, — XGMI/c is the spin of the MBH, 
with < < 1 the dimensionless spin parameter. The spin 
of the MBH evolves correspondingly by (King et al. 2005) 



dj./dt=-27r/ T^ivMR. 



(15) 



Since the torque is perpendicular to the MBH spin, it changes 
only its direction, not its magnitude. The BP torques work 
to align L and J,. The frame dragging torques dominate 
over the viscous torques out to the distance where the warp 
diffusion rate becomes faster than the Lense-Thirring pre- 
cession rate, /?bp = ^ V2 (/?bp ) / i^LT (^bp ) ■ Typically /?bp > 

Ra- The integrand in Eq. (fTsT i scales as T.{R)R 
so the torque on the MBH is dominated by contributions 
from the innermost region of the disk that remains non- 
aligned, and can be estimated by Tbp ^ ^^bp^bp(^bp) ^ 
?rV2(i?Bp)^(^Bp)VGM.i?BP- The timescale for the MBH spin 
to afign itself with the disk is then f|| ^ /.//bp, which in 
steady-state (M ~ 37rviE, V2 — Vi/2a^) can also expressed 
as 



-6xaUM,/M)Jrg/RjiP 



(16) 



The torqu e on the MBH spin by direct accretion of matter, 
Jncc ^ M\/GM,Ra, can be neglected relative to that by the BP 
effect, /bp. since /accZ-^BP ^ ^("-l \/ Ra / Rbv ^ 1- Accretion 
can no longer be ignored on timescales long enough for it to 
change the MBH mass appreciably, t>tE = rjcGr /4-nGmp = 
774.5 X 10^ yr, where is the e-folding, or Salpeter, timescale 
for growth by Eddington-limited accretion, Gt is the Thom- 
son cross-section, nip the proton mass and 77 the radiative ef- 
ficiency of the accretion. 

2.2. The internal structure of the disk 

The internal structure of the disk needs to be specified 
to determine the viscosity. In the a-disk prescription, the 
mid-plane temperature T reflects the azimuthal kinetic vis- 
cosity of the disk via Vi = a\c^ /Q-k, where cj = kT / ^ is 
the isothermal sound speed, pL is the mean molecular weight 

' The divergence of the BP torque toward the center is softened by substi- 
tuting R max{R, 300rj). 



and D.K = ^ GM, / R^ is the Keplerian angular frequency. The 
sound speed, in turn, determines the disk's exponential scale 
height H = Ci/Q.K. 

The physical parameters of the masing region in the ac- 
cretion disk of NGC4258 can be inferred from the observed 
aspect ratio H /R ^ 0.002, wh ich suggests a temperature of 
T - 600K d Argon et al.1 120071) for a gas pressure supported 
dis k (obs ervations disfavor magnetic support, see review by 
[Lo1|2005[) , and from the fact that H2O maser emission is pos- 
sible there. This requires the mid-plane molecular hydrogen 
density «Ht to be in the range lO^cm-3 to 10" cm (all the 
hydrogen is assumed to be molecular), the mid-plane gas tem- 
perature T to be in the range 300- 400 K to lOOOK, and 
the mid-plane pressure p/k to be in the range lO'^Kcm^^ 
to lO'^Kcm"^ ( lMalonevll2002l) . In addition, stabiHty against 
fragmentation requir es that Toomre 's criterion holds localljQ, 
Q = CiQ.K I TtGl, > 1 (iToomrdl 19641) . which further constrains 
the disk mass and temperature. 

It is unclear whether the accretion disk's internal viscosity 
alone can provide eno ugh heat to warm the molecular gas to 
masing temperatures. iNeufeld & MalonevI ( [T995ft argue that 
heating by X-ray irradiation by the central source is essential 
for creating the required conditions in the masing region. Out- 
side the masing region, the disk's properties are not usefully 
constrained by current observations. 

Here we make a choice of convenience to adopt the known 
solutiorQ of a gas pressure-dominated thin accretion disk, 
where the temperature is completely determined by the in- 
ternal viscosity (that is , no external heating) an d free-free 
(Kramer's law) opacity (IShakura & Sunvaevll973h . The pref- 
actor Ka of Kramer's law, K = Ka{p / Pc^{T /Ta)^'' ^'^ , where 
p = Y./ {\/2nH) is the total mid-plane mass density (the sub- 
script a denotes values at Ra) is then adjusted to be high 
enough to maintain the required temperature at the masing 
regiorQ. We assume Solar abundances {X — 0.7057, lArnettI 
119961) and that all the hydrogen is molecular, which corre- 
sponds to a mean molecular weight of jx = 2.358mp. The 
adopted value [H /R]a = 0.002 fixes the surface density E„ = 
PaV2n[H/R]aRa, where p„ = nii^_{Ra)2mp/X. We find that 
in order to have masing conditions between Ri, and Ri,, it is 
necessary to assume the highest temperature value at the in- 
ner edge, Ta ~ 10^ K. The assumed density there, «//2(^a) ~ 
3 X 10^ cm^-', which is within the masing range and consistent 
with the limits on pa/k — [pa/ li)Ta, is chosen so the total disk 
mass is 3 X IO^Mq. The mid-plane temperature is given by 
GsqT^ = f27/32)E^n|vi K (e.g. IP^aiik et al.l2002l Eg. 5.76), 
which here yields 



T{R,t) 



27k 



1/2 



aiKaTc 



7/2 



-I.'{R,t)D.liR) 



1/7 



(17) 



/27r32(7sB M'/^ P" 

For the values of Ta and Pa above, Eq. ( flTl i requires for 
self-consistency (Qi = 7.195 x 10^/ai cm^g^'. This is an 

^ Toomre's criterion can be recast as an approximate global condition, 
M(< R)/M. < a/vK = H/R, where vk = y'GM./R. 

^ The steady state solution has T « R^^/", £ « ^. ^ ^ ^ 

R'/*, oc i?-lV8_ vi 0= V2 « R+^''^, p/k oc R-21/8 and Toomre's Q « R-'/'* 

^ The actual opacity law in a co ol molecular disk presu mably resembles 
that of a proto-planetary disk, e.g. ISemenov et alj <2003l) . which does not 
have a simple analytic form. 
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extremely high near/mid IR opacity (Amax ~ 3 — 10 jim for 
T ^ 300-1000 K), as compared, for example, to the mean the- 
oretical opacity of (K-(2.2jUm)) ^ 380 0±700cm^gy ^ used to 
model dusty star forming cores (Sh irley et al.ll201 ll) . This is 
in line with the conclusion that an external source of heating 
is required to explain the gas temperature. 

2.3. Numerical implementation 

The finite difference method is used to integrate the evolu- 
tion equations (Eq. [TT]with G = c = M, = \) for the angular 

momentum surface densities {L; at time tj, at positions 

^ logarithmically spaced grid between R\ = Ra at 
the innermost stable circular orbit (ISCO) and Rn =Rp ^ R\- 
The edge points Rq and R^+i are used to enforce the bound- 
ary conditions. A brief desc r iption of the n umerical scheme 
jPapaloizou & Pringld [T983I: iPringld [T992b IS given m ap- 
pendix [A] 

2.3.1. Initial conditions and boundary conditions 

The initial conditions are a flat disk with surface density 
E = 'La{R/Ra)^^^'^ and normal £q, where E is chosen so as to 
satisfy the maser conditions between the masing region limits 
Ra and R^. The MBH spin is initially along the z-axis. We 
assume the boundary conditions L{Ri) = and d£/dR\R^ = 
(central mass sink and no warp), and d{viL,) /dR\R^ = (no 
torque). These boundary conditions allow mass through the 
inner boundary into the MBH. To prevent the disk from drain- 
ing on the viscous timescale, and from drifting away from the 
masing conditions much before that, a source term is added at 
the outer edge. It is adjusted every time-step tj to compensate 

for the mass lost from the disk, {AMy = IJLi(m/ -M/"'), 

where M/ = 27zRiARiT.j is the mass in a disk ring of width 
ARi = Ri — Ri-i <C Ri- The angular momentum density at 
the outermost grid point is then augmented by an amount 



{ALNyto, such that 



E-' - 



(AEa 



/Ra 



where {AL^y = —XMiAMy/lnRNARN (Xm is an order unity 
stabilization factor, see below). The magnitude of the angular 
momentum needed to compensate for the mass loss is there- 
fore 

(AL^)^' = -4-4 + 



[4- 



-(AE^).']2/?^-(L^)2.(18) 



The fluctuating factor Xm = 1 ± £ stabilizes the disk against 
the accumulation of numeric integration errors by allowing 
small over- or under-corrections, as needecQ. When {AMy < 
0, Xm is set to 1 + e, and conversely, when (AM)-' > 0, it is 
set to 1 — e. Experimentation indicates that e = 0.1 is a good 
choice. 

We find that with the source term thus defined, a flat 
disk with the inner boundary conditions ViEj/jj = 0, rapidly 
converges to a close approximation of its theoretically ex- 
pected steady state solution with the expected mass-loss rate 
of M ~ 37rviE (for R ^ Ra). Likewise, the mass of a non- 
stationary disk rapidly converges to a steady state value that 
is typically within (9(10^^) of its initial mass, even as its 
geometry continuously changes. The mean accretion rate 



over some period fy, to tj2 can then be estimated by ^M) = 

(If=,- Xm(AM)>)/(02-0i). 

This mass replenishment scheme represents an idealized 
case of a self -regulating mass supply, which prevents both the 
draining and over-loading and fragmentation of the disk. It 
should be emphasized that this assumption is introduced here 
for convenience only, to stabilize the disk long enough to col- 
lect robust statistics. Our statistical conclusions about the dy- 
namical response of the disk to the external stellar torques do 
not depend strongly on the existence of such a mass source, 
as long as the disk mass is in quasi-equilibrium for at least the 
relatively short RR timescale. 

2.3.2. The stellar resonant relaxation torques 

The RR torques due to the O(IO^) stars in the radius of 
influence cannot be modeled directly. Instead, they are ap- 
proximated here by representing the stellar cusp by a small 
number of concentric spherical shells delimited by the radii 
{>'k}^Lo (where ro = 0). The spacing between consecu- 
tive shell radii is chosen to be large enough so that the 
residual forces are approximately independent of each other, 
fk+i/fk > 22/(^^1'), which corresponds to the requirement 
that ^A^*(< rj,) > 2y^A?*(< r^_i) (BA 09), and is broadly 
consistent with the rigorous derivation of lKocsis & Tremaind 
(1201 ll) . The maximal number of such logarithmically spaced 
independent star shells in the simulations of NGC 4258 that 
are presented below is typically Ns = 5. The vector RR 
torques from each shell are represented by the gravitational 
field of a thin ring of mass M/^ = -yA^^^^T^J^^A^i^^T^TjM* 
with radius r^, = (rj^ + ri^^i)/! and normal n/.(f), which is 
interpolated smoothly in time by cubic spline from a se- 
quence of isotropic random normals {n^}^'™^''''''^' at times 
tj = jto{rii), where fjirn is the duration of the simulation and 
fo(r) = min(fsq(r),freact('")) (Eqs. |2]|9l) is the coherence time. 
This simulates the random reorientation of the local residual 
force after a vector RR timescale. Finally, the torque exerted 
by a star ring at r^, on a disk gas ring at which is directed 
along the line of intersection of the two rings, is calculated nu- 
merically by integrating over the gravitational force between 
all mass elements in the rings (iNayakshini,2005T) . 



MiMkRifk sinj3 

r^"" sin^isin02 
Jo [1 - dcosAj^/" 



(19) 

(£,xn,),(20) 



where M,- is the mass of the disk annulus, /3 = cos ' (£,■ • n^.), 
5 = 1- {rj - Rj) / {rl + Rj) anc0cosA =cosj3sin0i sin ^2 + 

COS^l COS 02- 

2.3.3. Validation of the code 

We verified that the numeric integration scheme conserves 
angular momentum to floating point precision, and conserves 
mass to a fractional precision of 0{lQ-^°yr-^), by switch- 
ing off the source term and evolving a flat disk from an out- 
of-steady-state initial surface mass density profile (E °^ R^^). 
We also evolved an E R^^^'^ initial configuration with a 



^ This algorithm of stabilization by over-shooting is inspired by the engi- 
neering method of control system hysteresis, used e.g. in thermostats. 



* The divergence of the torque when — > R, is smoothed by substituting 

5 ^ 5 = 1 - max [{rl -R}),{h - h-O^iR, - R,-lf] /(rf + 
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source term at R2 and with boundary conditions ViE = 
at and verified it approximates the analytical solution 
37rViE = M(l - ^R\/R) (e.g. iFrank et al.|[2002l) reasonably 
well (the match with the analytic solution is better on a linear 
grid than on a logarithmic one; overall the match o n a log- 
arithm ic grid is at a level similar to that obtained by iPringld 
II992I) . We also reproduce qualitatively the response of the 
di sk to initial twists and to the BP torques that are presented 
bv lMngle(199Z) . 

3. RESULTS 

We carried out a suite of disk evolution simulations, us- 
ing a stellar cusp model based on NGC4258 (Section 1 1.1b 
on a logarithmic grid with with = 100 points extending 
from Ri = 6rg to R^ = 1.5 x lOV^, (r^ ~ 5.5 x lO'^cm ~ 
1.8 X lO^^pc). The typical time-step increments in our sim- 
ulations were in the range Atj ^ 0.01-1 yr The initial disk 

configuration was flat with a E /?^^/^ surface mass density 
profile, and the disk was tilted by an angle /q relative to the 
z-axis, which coincided with the initial MBH spin, in those 
simulations with X ¥'^- I" the simulations presented below, 
we explore the initial conditions ;o =0,3/4 rad and x = 0, 1. 
At later times, we generalize the definition of the tilt angle 
to = cos"' (J, • Jdisk/-^«-^disk). where Jjisk is the total angu- 
lar momentum of the disk. To isolate the effect of the BP 
torques on the disk, we ran some simulations with A^* = 0. In 
the presence of RR, the disk never reaches a steady state, and 
the simulations are terminated after they are observed to reach 
statistical stationarity (recall that the mass supply rate is ad- 
justed continuously to maintain constant total mass. Section 
12.3.1b . Simulations without RR are stopped when the relative 
rate of change in {L,}^j falls below some very small value. 

Figure shows snapshots from a simulation of a mis- 
aligned disk (/() = 3/4 rad) evolving under both BP and RR 
torques, and the evolution of the inclination angle / for this 
and two other simulations, one with no RR and the other start- 
ing with an aligned disk (/q = rad). The evolution from 
the initial (and artificial) misaligned configuration, where the 
disk is tilted relative to the MBH spin all the way down to 
shows strong stochastic RR-dominated evolution that is su- 
perposed on a slower frame dragging-dominated secular evo- 
lution, until at later times the stochastic behavior dominates. 

3.1. Disk warping 

Figure ^ shows that RR-induced warping can easily re- 
produce the o bserved 8° warp angle across the maser region 
of NGC4258 (iHerrnstein et al.ll2005T) . without requiring any 
particular initial conditions. The three scenarios explored 
here, RR with maximal (X = 1) frame-dragging and a large 
(;o = 3/4 rad) initial tilt, or with a small (;o = rad) initial 
tilt, or without frame-dragging (x = 0), all exhibit a proba- 
bility Pa ^ 0{l) to have o > 8° (see probabilities quoted in 
Figure |2l). This should be contras ted with the scenario where 
only frame dragging is operating jCaproni et al.ll2007l) . In the 
example shown here, the warp angle never exceeds 5°. By 
careful choice of the initial parameters, it may be possible to 
obtain a large enough warp angle, and satisfy the constraints 
set by the direction of the radio jet with frame dragging only. 
However, this requires considerable fine-tuning, and may re- 
quire a very long-lived disk for the warp to grow ( iMartinI 
l2008h . 

The warping probability increases with the steepness of the 
stellar cusp, as more stars are available near the disk to torque 



it. For example, Pg = 0.02, 0. 1 1 and 0.29 for /q = rad = 1 
in a sequence of models with 7=3/2, 7/4 and 2, respectively. 

We find that for the /q = rad models of NGC4258, the 
rms RR-induced misalignment between the MBH spin axis 
and the disk's total angular momentum is rms(/) = 15° after 
10^ yr for the case X ~ as large as rms(;) = 44° for the 

X = case. This is an example of the stabilizing effect of the 
BP torques, which introduce a preferred plane for the disk, 
and suppress warps. It is interesting to note that the observed 
maser disk in NGC4258 is tilted by 30° to the present jet 
axis, as identified by the N orth and South hotpots dCecil et al.l 
120001: lWilsonetal.1120011) . Assuming that the jet is ahgned 
with the MBH spin axis, then such a tilt is consistent with 
RR torquing of a disk (see Figure[T]i), possibly around a sub- 
maximal spinning MBH. 

The warping of the disk exposes it to ionizing radia- 
tion from its innermost parts, which may play a central 
role in determining the physical conditions in the outer 
regions of the disk. The fraction of the central lumi- 
nosity that falls on the disk (assuming isotropic emis- 
sion) is expressed by the disk's covering fraction, Cf = 

/o'^d^ |cOS0niax(0) -COS0min(0)| /47r, where 0min and 0niax 

are the minimal and maximal inclinations above and below 
the disk's mean plane along R in azimuthal direction <j). Fig- 
ure ([3]) shows that RR-induced warping gives the disk a min- 
imal varying covering factor of C/ ~ 0. 1 — 0.3, on top of any 
contribution from a pre-existing large scale warp (e.g by the 
BP effect). 

3.2. Mass accretion rate 

As noted by iLodato & Pringld (l2006h . warping increases 
the mass accretion rate over that in a flat disk because of the 
additional dissipation due to the vertical shear We confirm 
here that the warps are associated with mass accretion fluc- 
tuations of up to factors of 3-4 over that in a flat disk (Fig- 
ure |2]i, resulting in an overall increase by a factor of up to 
ff^^lA in the average mass accretion rate. The power spec- 
trum of the accretion rate fluctuations in the maser model is 
very broad (Figure IH, and can be roughly approximated by 
a broken power-law, which flattens beyond > 10^ yr. This 
timescale corresponds to the RR coherence time in the outer 
half of the disk, r < R^ = 0.16 pc to r - /?2 = 0.27 pc. A 
rough estimate of this variability timescale can be obtained 
by evaluating the coherence time at the disk's mass-weighted 
mean radius, since that is where the gravitational coupling 
with the stars is maximal. That length-scale is typically in 
the regime where the back-reaction time is shorter than the 
self-quenching time (Eq|9]l, and so 

fvar freact(/?./) ^ {M,/Md)P{Rd) ~ 1.1 X 10^ (21) 

A similar correspondence is seen between the variability 
timescale of the AGN model (Section [3.3b . fvar — 1240, and 
the power-law break in its power spectrum (Figure I?]). Thus, 
fvar can be interpreted as the shortest timescale for which there 
is substantial variability power. 

3.3. MBH spin evolution 

The enhanced mass accretion rate translates directly to a 
faster rat e of chan ge in the magnitude of the the MBH spin, 
y, ~ My/GMtRa, as the MBH spin vector remains aligned 
with the inner regions of the disk due to the BP torques. The 
BP torques also couple the fluctuating RR torques on the disk 
to the MBH spin orientation (Figure |5] [TSl ). For the effect 
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Fig. 1. — Top row and bottom left: Snapshots from a disk evolution simulation of NGC4258 with ai = 0.25, X = Kramer's opacity law, and the initial 
conditions J, = z and ('o = 3/4 rad. The arrow denotes the MBH spin vector. The masing region is marked in blue. Top: The intersection of the disk mid plane 
with the X — z plane (the simulated plane of sky). Bottom left: The disk mid plane in 3D. Bottom right: The evolution of cos;' in the simulation shown in the 
snapshots, and for similar simulations with i'o = rad, and with BP torques only (W^, = 0). 

"catch up" with the faster changes in the direction of the stel- 
lar torques. We find that for this specific model, the MBH spin 
scatters across ~ 5° over an e-folding time (4.5 x 10^ years). 
Another consequence of the short coherence time (compared 



to be significant, the disk must be massive enough to carry a 
substantial amount of angular momentum (Eqs. [14] [TsT i, but 
not so massive as to remain effectively immobile against the 
RR torques. 

We find that the response of the MBH spin in NGC 4258 to 
the BP torques is very small, because of the very low mass of 
the disk. However, the effect can be larger in systems with a 
more massive accretion disk. Figure © shows as an exam- 
ple the evolution of the MBH spin axis relative to its initial 
direction (we omit the much smaller spin evolution due to ac- 
cretion), and that of the total angular momentum of the disk, 
for a low-mass luminous AGN model that corresponds to a 
Seyfert galaxy (L ~ 7 x lO'^^ergs^' for a radiative efficiency 
of T] = 0.1). The angular momenta of the MBH and the disk 
both execute a random walk around their initial orientation. 
The amplitude of the disk angular momentum jitter is larger 
than that of the MBH spin, because the RR coherence time, 
freact = 620 yr at R^/ = 2.2 x lO^^* pc is much shorter than the 
BP timescale fii = 4.1 x 10^ yr, and so the MBH spin cannot 



to fwaip(^rf) — 3.5 X lO"* yr) is that the RR-induced warps in 
this disk model have angles of only ~ 1°. Their effect on the 
mass accretion rate is correspondingly small, only a ^ 2% in- 
crease over that in a flat disk. 

4. DISCUSSION AND SUMMARY 

4.1. Discussion 

A circumnuclear accretion disk does not exist in isolation 
around a MBH. Rather, it shares the volume with the high- 
density nuclear cluster that is predicted to form there while the 
nucleus evolves. The cluster is expected to be on average spa- 
tially isotropic near the MBH, but the Poisson fluctuations in 
the distribution of orbital inclinations lead to a slowly varying 
residual force that exerts coherent torques on the disk (RR). 
We argue (Section [1.31 ) that a thin Keplerian disk is primar- 



Bregman & Alexander 



5-10 
5 4-10" 

o 

^ 3-10 

T3 



-5 



-5 



2-10" 



1-10 



-5 



— 1 1 1 1 1 1 1 1 1- 

BP only (io=3/4 rad) [<dM/dt>= 1.61e-05 Mp/yr, x1.19] 
RR only [<dM/dt>= 1 .59e-05 Mg/yr, x1.17] 
BP+RR (io=3/4 rad) [<dM/dt>= 1 .90e-05 Mg/yr, x1.40] 
BP+RR {\q=0 rad) [<dM/dt>= 1 .68e-05 Mg/yr, x1.24] 

Flat disk 




20 



D5 
CD 

O) 

C 

CO 

Q. 



1 1 1 1 1 

BP only (lo=3/4 rad) [Pg^O.OO] 

RRonly [Pg^O.IQ] 

BP+RR (io=3/4 rad) [Pg=0A4] 

BP+RR (io=0 rad) [Ps^O.H] 

NGC4258 angle 




BP only (io=3/4 rad) 
RR only 
BP+RR (if|=3/4 rad) 
BP+RR (io=0 rad) 



4-10' 



Fig. 2. — Simulated disk evolution with ai = 0.25, = or 1 and with Kramer's opacity law, for the cases of BP warping only {N^, = 0), RR warping only 
(X = 0) both BP and RR warping (all models have the same RR realization and initial MBH spin J, = z and iq = or 3/4 rad. They were evolved over lO' yr, 
of which the initial 4 x 10** yr are shown). Top: The mass accretion rate. The mass loss rate from the equivalent steady state flat disk, M = 1.4x lO^^M^yr^', is 
also shown for reference. The mean mass loss rates, and their enhancement over that of a steady state fi at disk are indicated in the plot labels. Bottom: The warp 
angle across the maser zone. The best fit warp angle of Q) = 8° for NGC4258 iHermstein et al.ll200^ is shown for reference. The fraction of time the models 
spend with a warp as large as that observed, Pg, are indicated in the plot labels. 

(Figs. [U 111, its mass accretion rate (Figs. |2l ^ and cover- 
ing factor (Figure (3]). The strong coupling between the stellar 
potential fluctuations and the disk via RR is then further ex- 
tended to the MBH spin via the frame-dragging BP effect, 
and this allows angular momentum to be transferred from the 
nuclear cluster to the MBH. The combined effects of the per- 
pendicular RR and BP torques excite a jitter in the MBH spin 
direction (Figure |5]l. In addition, the increased mass accre- 
tion rate due to the RR-induced warping, together with the 
disk / MBH spin alignment due to the BP effect, lead to an 
increased growth rate of the MBH spin magnitude. We con- 
clude that gravitational interactions between the stars and the 
disk excite a substantial level of irreducible variability in the 
disk properties — stationary accretion in a circumnuclear disk 
is merely an idealization. 

The large-scale warping of the disk (for NGC 4258, this 
coincides with the maser region warp) reflects the RR co- 
herence timescale at the mass-weighted mean disk radius, 
fvar {M,/Md)P{Rd) yr (- 10^ yr for NGC 4258). Gener- 
ally, the effect of RR on the disk and the MBH will be sub- 
stantial when the disk is long-lived, fdisk > fvar- There is a 
large spread in the estimates of AGN and QSO lifespans in the 
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Fig. 3. — The evolution of the maser disk's covering factor due to BP and 
RR for the models shown in Figure 0. 

ily torqued by these purely gravitational interactions with the 
stars, rather than by the hydrodynamical ones that occur when 
stars plunge through the disk. The RR torques warp the disk 
and can lead to order unity variability in the disk geometry 
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Fig. 4. — Left: The smoothed Normahzed Lomb-Scargle Periodograms (power spectra) of M(t), as function of inverse frequency, for the NGC 4 258 m aser 
disk model and for the AGN model discussed in Section f33l The coherence times of the star shells that are used to simulate the RR torques (Section [2.3.2) for 
the maser disk model are marked by crosses. The periodograms were smoothed to highlight the power-law break at ty^ (green line and circle), which indicates 
that most of the power is at ; > fyai = (M,/Mj)P(Rj) (The AGN model is displayed shifted on the logarithmic scale for comparison with the maser disk model). 
Right: The correlation between the variability on the large (maser region) scale and the relative mass accretion rate (normalized to the steady state case) in a 
model of NGC 4258 (BP and RR, ('o = rad and % = !)■ The mass accretion rate curve was smoothed by cubic spline, to filter out the short time-scale variability 
(cf Figure|2) and highhght the correlation between warping and mass accretion rate on longer time scales. 
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Fig. 5. — The evolution of the MBH spin J., and the disk's total angular 
momentum Jjisit. over an e-folding time (for 7] = 0. 1), shown projected on the 
(,v, y) plane, under the influence of BP and RR together, for a low mass AGN 
model with M. = 4 x 10'' Mq, = 2.3 pc, = 2, 7 = 2, ai = 0.1, X=l, 
Kramer's opacity law, normalized with \H/R\a = 0.002, p„ = 2.4 x 10"'^, 
T„ = lOOOK and = 10cm" g^' at = 0.1 pc, and the initial conditions 
J, =z and i = rad. These parameters result in = 1.6 x M,:,, Rj = 
2.2 X 10"' pc, ?var ~ 1240 yr, M = 0.013M,; yr"' = iArjME, Rbp = 3.6 x 
10"'* pc, and ?|| = 4. 1 x 10^ yr. The blue and red circles denote the projection 
of J, and Jdisk, respectively, at the end of the simulation. 



luminous p hase, fewx 10^-10^ yr (e.g. iGrazian et al.ll2004t 
iPorciani et al. 2004), with some recent analyses suggesting an 
even l onger lifespan of 5. 10^ yr (iGilli et al.ll2009l: iRoss et aTl 
120091; [Schawinski et al.ll2009l) . For long-lived disk systems, 
the cumulative effects of RR-torquing can be large. 

One implication of the factor ^ 1 .2 — 1 .4 increase in the 
mean m ass a ccretion rate over that in a flat stationary disk 
(Section [X2] l. is that an MBH fed by an RR-torqued disk, 
whose mass accretion rate is raised to the Eddington limit by 
warping, will be exp[(/^ — l)(rdisk/f£)] more massive than 
one that is fed by an otherwise identical flat disk, which is 
accreting at only l/fj^ of the Eddington Umit. For exam- 



ple, RR-torquing can accelerate MBH growth over a time 
?disk = 10f£ = 4.5 X 10** yr (assuming 77 = 0.1) by a factor of 
up to ^ 50 (for = 1 -4). Alternatively, if radiation pressure 
prevents the increased RR-induced mass flow from accreting 
on the MBH, the accumulating excess mass may trigger out- 
flows, or disk fragmentation and star formation. 

Another cumulative effect of RR-torquing of the disk is the 
displacement of the MBH spin from its initial orientation due 
to the BP-induced random jitter It is interesting to note that 
the few degrees amplitude of the displacement is still consis- 
tent with the typical obse rved ~ 5° opening angle of A GN 
radio jets on large scales dOppenheimer & BirettalfT994l) . as- 
suming the jet direction reflects the spin direction. 

A direct geometric consequence of RR-induced warping 
is the substantial time-averaged covering factor the disk ac- 
quires relative to the central continuum source. This allows 
the disk to intercept the central radiation and be heated and 
ionized by it. Such X-ray irradiation may in fact be essen- 
tial for raising the gas temperature to the ra nge necessary for 
maser emission (Neufeld & Malone"vl ll995h . for the produc- 
tion of the observed e mission lines in AGN disks in general 
(ICollii>S ouffrinlll9 87|), for dr iving winds and outflows from 
disks (jProga & KaII manl2004l) . or for explaining Narrow Line 
Seyfe rt 1 galaxies and ultra-soft AGN dPuchnarewicz & Sorial 
I2OOI . 

The magnitude of the effects of the stars on the disk de- 
pend on the central concentration of the stellar cluster, and 
is therefore sensitive to the degree of mass segregation in the 
system and to the fraction of massive compact remnants in 
the population. The more centrally concentrated the stellar 
density profile, and the more massive the stars (for RR, the 
relevant quantity is (M^) / (M*), iRauch & Tremaind [19961) . 
the larger are the effects on the disk. It is worth noting 
that strong mass segregation with a power-law cusp pro- 
files of 7 > 2 is expected to occur in nuclei with old stellar 
populations ([Ale xander & Ho pinarill2009l: fKeshet et al.ll2009l: 
iPreto & Amaro-S eoane 2010), where the inner parts of the 
cusp are dominated by stellar black holes of mass 0{IQM,7)). 
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It is thus plausible that a substantial fraction of galactic nu- 
clei have conditions that are conducive to RR-torquing of an 
accretion disk. 

In this study we presented numerical results for galactic 
nuclei with MBH masses in the range ^ 4 x 10^ Mp, (with 
a massive disk close to Me) to ^ 4 x lO^M© (with a low- 
mass, low-accretion rate disk). In order to scale the results 
to other systems, it is necessary to specify the scaling of the 
disk properties with MBH mass. We defer this to future work, 
and present here only a simple preliminary analysis that sug- 
gests that the effects are generally more important in lower- 
mass MBHs. Assuming that all AGN disks extend up to some 
universal large multiple of rg (e .g. ~ 2000ri, for a-disks lim- 
ited by gravitational instability, lGoodmanll2003h . and assum- 
ing they all have the same aspect ratio, then R^i oc M, and 
Mj/M, H/R, so that <^ M,. The angular momentum 
in the disk therefore scales as J^i cx Mj^M,Rd °^ M^. The 
M,/a relation indicates that the MBH radius of influence 

scales as r/, M,'"'^ (assuming M, a'*), so that the num- 
ber of stars in the MBH cusp inside R^i scales as Ni,{Rd) °^ 

Mt{Rd/rh) °^ M, . The residual angular momentum 
in the stars that are efficiently coupled to the disk then scales 

as Jn °= ^N^{Rd)VM^ °^ mI'^'^^^'^, and it then follows that 

Jj/Jn °^ The ratio of the disk and stellar angular 

momenta is a slowly rising function of M, for 7 > 1 . This 
implies that, all other parameters being equal, it should be- 
come progressively more difficult for the stars to torque the 
disk as the MBH mass increases. In addition, as argued in 
Section ( |1.3t , the larger M,, the more the disk responds to the 
RR torques as a rigid body, rather than by growing warps. 

4.2. Caveats 

The calculations and conclusions presented here are limited 
by various assumptions and approximations. They are listed 
here briefly. Some of these issues will be addressed in future 
work. 

Limitations of the galactic nucleus model — The RR torques by 
the stellar cluster are approximated by th e gravi tational field 
of a small number of thin rings (Section 12 3.21 ). which pro- 
vides only a rough approximation of the true power spectrum 
of the stellar perturbations. The back-reaction of the disk on 
the stars is approximated simplistically by limiting the RR co- 
herence time so it does not exceed the back-reaction time. A 
more realistic fluctuation spectrum will allow, among other is- 
sues, to better model the very short timescale fluctuations, and 
explore whether these have any relation to the observed opti- 
cal/UV variability of AGN. Our conclusions on mass accre- 
tion rates and MBH spin evolution are generalizations based 
on extrapolating simulations that were constructed explicitly 
to represent the disk and cluster of NGC4258. They should 
be scaled and tested for a wider range of MBHs, disks and 
nuclear clusters . 

Limitations of the disk model — The physics underlying ac- 
cretion disks are still not well understood. In particular 



the nature of the azimuthal and vertical visc osities, and the 
relation between them, are quite uncertain dOgilvid Il999t : 
iLodato&Pringld 120071: iLodato & Pric3 120101) . The disk is 
not modeled self-consistently. The thermal structure of the 
disk is based for simplicity on a Kramer opacity law with a 
free normalization constant, which is clearly non-physical, as 
indicated by the very high opacity that is required to produce 
the masing conditions. External X-ray heating, which is likely 
important, is neglected. MBH spin evolution due to mass ac- 
cretion can be relevant on longer timescales, but is not taken 
into account, and the artificial BP torque suppression on small 
scales that is needed to stabilize the results numerically, may 
both lead to an underestimation of the effect of disk warping 
on the MBH spin evolution. Finally, the disk is approximated 
as Newtonian and Keplerian, even near the ISCO, which is 
held fixed at its Schwarzschild value. 

Limitation s of the numerical scheme — The numerical scheme 
used here jPapaloizou & Pringlel[T983l:IPringldll992h can not 
describe azimuthal modes and as implemented, is limited to 
moderate warp angles (although visco sities for arbitr ary large 
angles can be calculated numericallv. lOgilvi^l 19991) . In par- 
ticular, disk flipping from co-rotation to counter-rotation, or 
disk destruction by large deformations, cannot be modeled re- 
liably. 

4.3. Summary 

We analyzed and simulated numerically the evolution of a 
thin accretion disk around a MBH that is surrounded by a stel- 
lar cluster We took into account the disk's internal viscous 
torques, the frame-dragging torques of a spinning MBH and 
the stellar orbit-averaged gravitational torques. We show that 
the evolution of the MBH mass accretion rate, the MBH spin 
growth rate, and the covering fraction of the disk relative to 
the central ionizing continuum source, are all strongly cou- 
pled to the stochastic fluctuations of the stellar potential via 
the warps that the stellar torques excite in the disk. These 
lead to fluctuations by factors of up to a few in these quanti- 
ties over a wide range of timescales, with most of the power 
on timescales > {M, / M([)P{R^i). The response of the disk is 
stronger the lighter it is and the more centrally concentrated 
the stellar cusp. We demonstrated these effects by simulating 
the evolution of the maser disk in NGC 4258, and show that 
its observed (9(10°) warp can be driven by the stellar torques. 
We also show that the frame-dragging of a massive AGN disk 
couples the stochastic stellar torques to the MBH spin and can 
excite a jitter of a few degrees in its direction relative to that 
of the disk's outer regions. 
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APPENDIX 

A. THE DISCRETIZED DISK EVOLUTION EQUATION 

The integration of the evolution equation generally follows the scheme of lPringld ( 119921) . but some details of the implementation 
are different. It is presented here briefly for completeness. 
All quantities are expressed in a system of units where G — c — M, — 1 . Denoting by Az the equal spacing of the logarithmic 
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grid, the grid points are at Ri — /?ie'' for / = 0, . . . + 1, where points Rq, Rn+\ represent the boundaries and = 6 is at 
the ISCO. The angular momenta densities at time tj, V,\ {i~ I,... ,N), are advanced in time to tj+\ ~ tj + Atj by 



Ati 



+1 
j 



1,1-1^1-1 

•j 



-1, 



?/)-vi,_i,,Lti,,x(^/-£t,)]} 



(Al) 



(A2) 



where L^,-^i 



(l/+L^i)/2, -i 



ij+i 



, and v/ 



n.i.i+l 



V^j, j)/2 for n =2,3. The index k is defined as 



k = i—l when the advective velocity V^^^ ■ > 0, and k = i when V^^^ ■ < 0, where 



adv,( 



1 





—ii-l 


2 




2Az 





(A3) 



for / = 1,...N. For / = 0, ; — 1 in Eq. ( IA3b is substituted by 0, and for / = + 1 , / + 1 is substituted by + 1 . The viscosities 
(and all other thermodynamic properties) are then updated via Eq. ( fTTI ). 

The inner boundary conditions, L(/?i) = 0, dl/dR]^^ ~ are enforced by setting Lq = (this causes a slight deviation from 
the analytic solution, since Lj 0, but is not zero identically), and — t\. The outer boundary condition, d{y\V,) / dR = is 
enforced by setting L^^j = (v/^/Vj ;y^j)L-(/. The mass-loss from the inner edge of the disk is balanced (in the statistical sense) 

by adding a source term at the outer radius, Lji^ + {ALnY £q (see Section |2.3.1| and Eq. [TSl l. 

The time-step size is adjusted every time-step to 

Atj = min(Ari^, , Af4,Af4,Af^p,Af^R)/2 , (A4) 

where 



Af4 = min,(A/?,)V<,i+i for« = 1,2,3 , 

Af^p = min,L//|T^p ,.| , Af^j, = min,(M,/M.) / ^N,{<Ri)R. 



3/2 



(A5) 



and where Vi .,.,+i is defined similarly to V2.,.,+i. 

Angular momentum conservation is monitored by evaluating the change in the total angular momentum of the disk be- 
tween times ry_i and tj in two different ways. One, which corresponds to the divergence term in the continuity equa- 



tion, is AJg(jgg(Lg,Lj,L^,L-^^j,Lg^ ,Lj^ ,1^^^ ,L^4i) = — 2;rEv+'/?,zi/;,(L;f ^ ), which can be written as a function 
of the disk edges only, since all the inner grid points cancel out in the sum. The other, which corresponds to the time 
derivative term in the continuity equation, is the change of angular momentum in the bulk, AJ^^jj. = Jbuik ^ Jbuik' where 
Jbuik ~ ^nHf^iRiARihj . The fractional degree of non-conservation per time-step in the absence of source terms (Section [2. 3. 3l l 



is then 5J^ = 



edge 



■^buik- Analogous expressions are used to monitor mass conservation. 



B. SECOND ORDER COEFFICIENTS OF THE VISCOSITY PARAMETER EXPANSION 

The expansion of the dimensionless viscosity parameters (Eq. fT2] | was derived by lOgilvi^ (119991) . The second order coefficients 
are reproduced here for convenience in the notation of this paper, as functions of the in-plane shear viscosity parameter ai and 
the adiabatic index F only (F = 7/5 for diatomic gas). The dependence on a bulk viscosity parameter, assumed here to be zero, 
is not included. 



-2(l-17af + 21af) 



B2 = 2Re(e23), Bi=lm{Q23), 



(Bl) 
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where 



and 



a + bT + c[T-i8/3)au] 
4ai[3-r+(8/3)ai/](-ai+2/)3(ai+2/) ' ^ ' 

:(12 + 49a?-408af + 2af) + (115 + 109ar-41af + 4af)ai!, 
{lS-S7af- 24-af) + (87 ~ 6af)aii, 

: ( 1 2 - 36ai2 + 140af + 2af ) + (25 - 1 1 + 2 1 af ) «! / . (B3) 
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